library(plotly)
library(tidymodels)
## install BiocManager if not installed
#if (!requireNamespace("BiocManager", quietly = TRUE))
#install.packages("BiocManager")
## then install mixOmics
#BiocManager::install("mixOmics")
library(mixOmics)
df <- read.csv('data2_areaR.csv',header = T, sep = ";",dec = ".")
rownames(df)= df[,1]
data <- df[,4:64]
categoria <- df$categoria
amostras <- df$amostra
bandeira <- df$bandeira
prin_comp <- prcomp(data, center = T, scale. = T)
components <- prin_comp[["x"]]
components <- data.frame(components, categoria, amostras)
#plotPCA
font.config <- list(
family = "sans serif",
size = 12,
color = 'black')
fig <- plot_ly(components,
x = ~PC1,
y = ~PC2,
color = ~categoria,
colors = c('#BF382A', '#0C4B8E'),
type = 'scatter',
mode = 'markers',
text = ~amostras)|>
layout(xaxis = list(title = "PC1 29,2%"),
yaxis = list(title = "PC2 17,8%"),
title = "",
font = font.config)
fig
#orca(fig, file="PCA.png", scale = 3)
#orca(fig, file="PCA.svg")
tot_explained_variance_ratio <- summary(prin_comp)[["importance"]]
tot_explained_variance_ratio
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 4.218463 3.294242 2.836302 1.905387 1.652718 1.416158
## Proportion of Variance 0.291730 0.177900 0.131880 0.059520 0.044780 0.032880
## Cumulative Proportion 0.291730 0.469630 0.601510 0.661030 0.705800 0.738680
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 1.38593 1.243424 1.236495 1.095549 1.050944 0.9942217
## Proportion of Variance 0.03149 0.025350 0.025060 0.019680 0.018110 0.0162000
## Cumulative Proportion 0.77017 0.795520 0.820580 0.840260 0.858360 0.8745700
## PC13 PC14 PC15 PC16 PC17
## Standard deviation 0.9692859 0.9302106 0.8494444 0.814257 0.7787514
## Proportion of Variance 0.0154000 0.0141900 0.0118300 0.010870 0.0099400
## Cumulative Proportion 0.8899700 0.9041500 0.9159800 0.926850 0.9367900
## PC18 PC19 PC20 PC21 PC22
## Standard deviation 0.6944917 0.6296314 0.5886719 0.5500549 0.5179903
## Proportion of Variance 0.0079100 0.0065000 0.0056800 0.0049600 0.0044000
## Cumulative Proportion 0.9447000 0.9512000 0.9568800 0.9618400 0.9662400
## PC23 PC24 PC25 PC26 PC27
## Standard deviation 0.4718777 0.4624003 0.4404313 0.4373824 0.3995965
## Proportion of Variance 0.0036500 0.0035100 0.0031800 0.0031400 0.0026200
## Cumulative Proportion 0.9698900 0.9733900 0.9765700 0.9797100 0.9823300
## PC28 PC29 PC30 PC31 PC32
## Standard deviation 0.3620183 0.3287083 0.3105873 0.3031722 0.2879618
## Proportion of Variance 0.0021500 0.0017700 0.0015800 0.0015100 0.0013600
## Cumulative Proportion 0.9844800 0.9862500 0.9878300 0.9893400 0.9907000
## PC33 PC34 PC35 PC36 PC37
## Standard deviation 0.2678198 0.2398007 0.2296439 0.2085803 0.2007504
## Proportion of Variance 0.0011800 0.0009400 0.0008600 0.0007100 0.0006600
## Cumulative Proportion 0.9918700 0.9928100 0.9936800 0.9943900 0.9950500
## PC38 PC39 PC40 PC41 PC42
## Standard deviation 0.1930121 0.1781942 0.1746784 0.1704367 0.1633052
## Proportion of Variance 0.0006100 0.0005200 0.0005000 0.0004800 0.0004400
## Cumulative Proportion 0.9956600 0.9961800 0.9966800 0.9971600 0.9976000
## PC43 PC44 PC45 PC46 PC47
## Standard deviation 0.1526067 0.1441266 0.1319377 0.130614 0.1289309
## Proportion of Variance 0.0003800 0.0003400 0.0002900 0.000280 0.0002700
## Cumulative Proportion 0.9979800 0.9983200 0.9986000 0.998880 0.9991600
## PC48 PC49 PC50 PC51 PC52
## Standard deviation 0.1204568 0.1072623 0.08702224 0.0731475 0.07191255
## Proportion of Variance 0.0002400 0.0001900 0.00012000 0.0000900 0.00008000
## Cumulative Proportion 0.9993900 0.9995800 0.99971000 0.9998000 0.99988000
## PC53 PC54 PC55 PC56 PC57
## Standard deviation 0.06117772 0.0365507 0.0263285 0.02460594 0.02035427
## Proportion of Variance 0.00006000 0.0000200 0.0000100 0.00001000 0.00001000
## Cumulative Proportion 0.99994000 0.9999600 0.9999700 0.99998000 0.99999000
## PC58 PC59 PC60 PC61
## Standard deviation 0.01637325 0.01056493 0.009337684 0.008251142
## Proportion of Variance 0.00000000 0.00000000 0.000000000 0.000000000
## Cumulative Proportion 1.00000000 1.00000000 1.000000000 1.000000000
tit = 'PCA 3d'
fig1_3d <-plot_ly(components,
x = ~PC1,
y = ~PC2,
z = ~PC3,
color = ~categoria,
colors = c('#BF382A', '#0C4B8E'),
text = ~amostras,
type = "scatter3d",
mode = "markers") |>
layout(scene = list(xaxis = list(title = "PC1 29,2%"),
yaxis = list(title = "PC2 17,8%"),
zaxis = list(title= "PC3 13,2%")),
title = "PCA")
fig1_3d
# orca(fig1_3d, file="PCA3d.svg")
fig2 <- plot_ly(components,
x = ~PC1,
y = ~PC2,
color = ~bandeira,
colors = c('orange','black','#BF382A','green4','#0C4B8E'),
type = 'scatter',
mode = 'markers',
text = ~amostras)|>
layout(xaxis = list(title = "PC1 29,2%"),
yaxis = list(title = "PC2 17,8%"),
title = "",
font = font.config)
fig2
# orca(fig2, file="PCA_band.png", scale = 3)
plot_ly(components,
x = ~PC1,
y = ~PC2,
z = ~PC3,
color = ~bandeira,
colors = c('orange','black','#BF382A','green4','#0C4B8E'),
text = ~amostras,
type = "scatter3d",
mode = "markers")%>%
layout(scene = list(xaxis = list(title = "PC1 29,2%"),
yaxis = list(title = "PC2 17,8%"),
zaxis = list(title= "PC3 13,2%")),
title = "PCA Bandeira")
Aqui, foi calculado uma PLS-DA levando em consideração a categoria para identificar qual são os marcadores importantes que distinguem cada marca.
X <- data
Y <- categoria
#Calculando a PLS-DA
result.plsda <- plsda(X, Y, scale = TRUE, ncomp = 5)
# CV
set.seed(123)
vald_plsda <- perf(result.plsda, validation = "Mfold", folds = 5, progressBar = F, auc = T, nrepeat = 100)
plot(vald_plsda, col = color.mixo(1:3), sd = T,legend.position="horizontal")
# https://rdrr.io/cran/mixOmics/man/plot.perf.html para ajudar no plot
# Escolhendo a melhor VL para a discriminação dos grupos
vald_plsda$choice.ncomp
## max.dist centroids.dist mahalanobis.dist
## overall 4 4 4
## BER 4 4 4
# AUC
vald_plsda$auc$comp4
## AUC.mean AUC.sd
## 0.968573000 0.003692385
# Plots
plsda_varcomps <- result.plsda$variates$X[,1:3]
datapls <- data.frame(plsda_varcomps,categoria)
cores <- c('#BF382A', '#0C4B8E')
# Plot interativo 2d
plot_ly(data = datapls,
x = ~ comp1,
y = ~ comp2,
colors = cores,
type = "scatter",
mode = "markers",
color = ~ categoria,
size = 1)%>%
layout(xaxis = list(title = "Componente 1"),
yaxis = list(title = "Componente 2"),
title = "PLS-DA")
# 3d
plot_ly(data = datapls,
x = ~ comp1,
y = ~ comp2,
z = ~ comp3,
colors = cores,
type = "scatter3d",
mode = "markers",
color = ~ categoria) %>%
layout(scene = list(xaxis = list(title = "Componente 1"),
yaxis = list(title = "Componente 2"),
zaxis = list(title= "Componente 3")),
title = "PLS-DA 3d",
font = font.config)
# Calculo VIPScores PLS-DA
vipscores <- vip(result.plsda) #calculo dos VIPScores da plsda
vipscores.df <- data.frame(vipscores,colnames(X)) #df juntando vipscores + nome das variáveis
vip_cutoff = 1.1 #parâmetro
# Calculo dos VIPScores mais importantes para quarta componente
vipscores4comp <- vipscores.df[which(vipscores.df[,4] >= vip_cutoff),c(4,6)] #esse comando seleciona com o which quais variáveis são maiores que o parâmetro '=>vip_cutoff = 1'. o 'c(1,4) significa, o 1 a primeira componentes e o 4 a coluna com o nome da variável. Foi escolhido a 4 comp no comando 'vald_plsda$choice.ncomp'.
# Lista de Variáveis mais importantes
vipscores4comp
## comp4 colnames.X.
## VOC10 1.247602 VOC10
## VOC11 1.307620 VOC11
## VOC13 1.190616 VOC13
## VOC25 2.837054 VOC25
## VOC28 1.170878 VOC28
## VOC29 1.175418 VOC29
## VOC31 1.158365 VOC31
## VOC32 1.166913 VOC32
## VOC33 1.171270 VOC33
## VOC34 2.054230 VOC34
## VOC36 1.110257 VOC36
## VOC37 1.190100 VOC37
## VOC38 1.182578 VOC38
## VOC40 1.216729 VOC40
## VOC41 1.210116 VOC41
## VOC42 1.195901 VOC42
## VOC44 1.695066 VOC44
## VOC59 2.184759 VOC59
A PLS-DA apresentou acurácia de 0,96857 com desvio padrão de 0,00369 na quarta variável latente. Excelente modelo para predição de EHA/EHC.
Foram selecionados na 4 componente, 27 variáveis com o VIPScore >= 1.
#Retirando o grupo shell da base de dados
df2 <- df[df$bandeira!='EHA-S',]
X2 <- df2[,4:64]
categoria2 <- df2$categoria
Y2 <- categoria2
#PLS-DA sem shell categoria
result.plsda2 <- plsda(X2, Y2, scale = TRUE, ncomp=5)
# CV
set.seed(1234)
vald_plsda2 <- perf(result.plsda2, validation = "Mfold", folds = 5, progressBar = F, auc = T, nrepeat = 100)
plot(vald_plsda2, col = color.mixo(1:3), sd = T,legend.position="horizontal")
#Selecionando a VL mais importante
vald_plsda2$choice.ncomp
## max.dist centroids.dist mahalanobis.dist
## overall 3 3 4
## BER 3 3 4
#AUC da 3 comp
vald_plsda$auc$comp3
## AUC.mean AUC.sd
## 0.963487000 0.003705924
#Plot 3d da PLS-DA Sem shell
plsda2_varcomps <- result.plsda2$variates$X[,1:3]
datapls2 <- data.frame(plsda2_varcomps,categoria2)
plot_ly(data = datapls2,
x = ~ comp1,
y = ~ comp2,
z = ~ comp3,
colors = cores,
type = "scatter3d",
mode = "markers",
color = ~categoria2) %>%
layout(scene = list(xaxis = list(title = "Componente 1"),
yaxis = list(title = "Componente 2"),
zaxis = list(title= "Componente 3")),
title = "PLS-DA 3d",
font = font.config)
#Calculo VIPScores PLS-DA sem Shell
vipscores2 <- vip(result.plsda2) #calculo dos VIPScores da plsda
vipscores.df2 <- data.frame(vipscores2,colnames(X2)) #df juntando vipscores + nome das variáveis
vip_cutoff2 = 1.1
#Calculo dos VIPScores mais importantes para primeira componente
vipscores3comp <- vipscores.df2[which(vipscores.df2[,3] >= vip_cutoff2),c(3,6)]
#Esse comando seleciona com o which quais variáveis são maiores que o parâmetro '=>vip_cutoff = 1'. o 'c(1,4) significa, o 1 a primeira componentes e o 4 a coluna com o nome da variável.
vipscores3comp
## comp3 colnames.X2.
## VOC10 1.266392 VOC10
## VOC11 1.289764 VOC11
## VOC13 1.176945 VOC13
## VOC23 1.141260 VOC23
## VOC25 2.974478 VOC25
## VOC34 2.301792 VOC34
## VOC44 1.973896 VOC44
## VOC48 1.163187 VOC48
## VOC51 1.116619 VOC51
## VOC52 1.278763 VOC52
## VOC59 2.305147 VOC59
Foi escolhida a terceira componente para a PLS-DA EHA/EHC sem o grupo Shell e 11 compostos foram selecionados com VIPScores >1.1 para essa componente.
# heatmap com melhor as variáveis da PLS-DA sem o grupo Shell
df.heatmap <- df2[,vipscores3comp$colnames.X2.]
df.heatmap_scaled <- scale(df.heatmap,center = TRUE, scale = TRUE)
bandeira2 <- df2$bandeira
df.heatmap.cat <- data.frame(bandeira2,df.heatmap_scaled)
variables <- colnames(df.heatmap)
values <- as.matrix(df.heatmap.cat[,2:12])
# cores
vals <- unique(scales::rescale(c(volcano)))
o <- order(vals, decreasing = FALSE)
cols <- scales::col_numeric("reds", domain = NULL)(vals)
colz <- setNames(data.frame(vals[o], cols[o]), NULL)
# heatmap bandeira
fig3 <- plot_ly(x=variables,y=bandeira2, z = values, type = "heatmap",colorscale = colz) |> layout(font=font.config)
fig3
# heatmap com melhor as variáveis da PLS-DA sem o grupo Shell
df2.heatmap <- df[,vipscores4comp$colnames.X.]
df2.heatmap_scaled <- scale(df2.heatmap,center = TRUE, scale = TRUE)
df2.heatmap.cat <- data.frame(bandeira,df2.heatmap_scaled)
variables2 <- colnames(df2.heatmap)
values2 <- as.matrix(df2.heatmap.cat[,2:19])
# cores
vals <- unique(scales::rescale(c(volcano)))
o <- order(vals, decreasing = FALSE)
cols <- scales::col_numeric("reds", domain = NULL)(vals)
colz <- setNames(data.frame(vals[o], cols[o]), NULL)
# heatmap bandeira
fig5 <- plot_ly(x=variables2,
y=bandeira,
z = values2,
type = "heatmap",
colorscale = colz) |> layout(font=font.config)
fig5
#orca(fig4, file="heatmap.png", scale = 3)
df.marcadores <- data.frame(categoria,df[,c(28,32,61)])
# Supondo que df.marcadores já foi criado conforme seu código
df.marcadores$categoria <- ifelse(df.marcadores$categoria == "EHA", 1,
ifelse(df.marcadores$categoria == "EHC", -1, df.marcadores$categoria))
# Transformando a coluna 'categoria' de volta para numérico
df.marcadores$categoria <- as.numeric(as.character(df.marcadores$categoria))
# Verificando a transformação
str(df.marcadores$categoria)
## num [1:197] 1 1 1 1 1 1 1 1 1 1 ...
#Regressão
lm_model <- linear_reg() %>%
set_engine('lm') %>%
set_mode('regression') %>%
fit(categoria ~ VOC25 * VOC29 * VOC59, data = df.marcadores)
summary(lm_model)
## Length Class Mode
## lvl 0 -none- NULL
## spec 7 linear_reg list
## fit 12 lm list
## preproc 1 -none- list
## elapsed 1 -none- list
## censor_probs 0 -none- list
#Predição
df_predict <- lm_model %>% predict(df)
rownames(df_predict) <- rownames(df)
## Warning: Setting row names on a tibble is deprecated.
df_predict
## # A tibble: 197 × 1
## .pred
## * <dbl>
## 1 1.15
## 2 0.477
## 3 0.762
## 4 0.808
## 5 0.970
## 6 0.883
## 7 1.48
## 8 0.444
## 9 0.909
## 10 1.15
## # ℹ 187 more rows
summary(lm_model)
## Length Class Mode
## lvl 0 -none- NULL
## spec 7 linear_reg list
## fit 12 lm list
## preproc 1 -none- list
## elapsed 1 -none- list
## censor_probs 0 -none- list
#Plot
fig6 <- plot_ly(df, y = ~ lm_model$fit$fitted.values, x = ~categoria,
colors = categoria,
type = 'scatter',
alpha = 0.2,
mode = 'markers',
name = 'Tips',
text = rownames(df_predict)) %>%
layout(xaxis = list(title = "Categoria"),
yaxis = list(title = "Valor Predito"),
title = "")
# Cálculo da acurácia
# Sendo, o positivo a aditivação:
# VP <- Qtd de EHA dado como EHA (qtd EHA > 0)
# VN <- Qtd de EHC dado como EHC (qtd EHC < 0)
# FP <- Qtd de EHC dado como EHA (qtd EHC > 0)
# FN <- Qtd de EHA dado como EHC (qtd EHA < 0)
df_predict1 <- data.frame(df_predict$.pred,rownames(df_predict), categoria)
VP <- df_predict1[which(df_predict1$df_predict..pred>0 & df_predict1$categoria=="EHA"),]
nrow(VP)
## [1] 91
VN <- df_predict1[which(df_predict1$df_predict..pred<0 & df_predict1$categoria=="EHC"),]
nrow(VN)
## [1] 98
FP <- df_predict1[which(df_predict1$df_predict..pred>0 & df_predict1$categoria=="EHC"),]
nrow(FP)
## [1] 0
FN <- df_predict1[which(df_predict1$df_predict..pred<0 & df_predict1$categoria=="EHA"),]
nrow(FN)
## [1] 8
AUC = 100*(nrow(VP)+nrow(VN))/(nrow(VP)+nrow(VN)+nrow(FP)+nrow(FN))
AUC
## [1] 95.93909
ESP=100*(nrow(VN)/((nrow(VN)+nrow(FP))))
ESP
## [1] 100
SENS=100*(nrow(VP)/((nrow(VP)+nrow(FN))))
SENS
## [1] 91.91919
#orca(fig6, file="fitted_values.png", scale = 3)
Das 8 amostras EHA que não foram classificadas corretamente, 4 eram de bandeira branca. Para testar vou adicionar um marcador para as “outras bandeiras”
df.marcadores2 <- data.frame(categoria,df[,c(28,32,61,37)])
df.marcadores2$categoria <- ifelse(df.marcadores2$categoria == "EHA", 1,
ifelse(df.marcadores2$categoria == "EHC", -1, df.marcadores2$categoria))
# Transformando a coluna 'categoria' de volta para numérico
df.marcadores2$categoria <- as.numeric(as.character(df.marcadores2$categoria))
# Verificando a transformação
str(df.marcadores2$categoria)
## num [1:197] 1 1 1 1 1 1 1 1 1 1 ...
#Regressão
lm_model2 <- linear_reg() %>%
set_engine('lm') %>%
set_mode('regression') %>%
fit(categoria ~ VOC25 * VOC29 * VOC59 * VOC34, data = df.marcadores2)
summary(lm_model2)
## Length Class Mode
## lvl 0 -none- NULL
## spec 7 linear_reg list
## fit 12 lm list
## preproc 1 -none- list
## elapsed 1 -none- list
## censor_probs 0 -none- list
#Predição
df_predict2 <- lm_model2 %>% predict(df)
rownames(df_predict2) <- rownames(df)
## Warning: Setting row names on a tibble is deprecated.
df_predict2
## # A tibble: 197 × 1
## .pred
## * <dbl>
## 1 1.02
## 2 0.526
## 3 0.963
## 4 0.794
## 5 1.02
## 6 0.878
## 7 1.35
## 8 0.543
## 9 1.05
## 10 1.11
## # ℹ 187 more rows
summary(lm_model2)
## Length Class Mode
## lvl 0 -none- NULL
## spec 7 linear_reg list
## fit 12 lm list
## preproc 1 -none- list
## elapsed 1 -none- list
## censor_probs 0 -none- list
#Plot
plot_ly(df, y = ~ lm_model2$fit$fitted.values, x = ~categoria,
colors = categoria,
type = 'scatter',
alpha = 0.2,
mode = 'markers',
name = 'Tips',
text = rownames(df_predict2)) %>%
layout(xaxis = list(title = "Categoria"),
yaxis = list(title = "Valor Predito"),
title = "PLS-DA")
# Cálculo da acurácia
# Sendo, o positivo a aditivação:
# VP <- Qtd de EHA dado como EHA (qtd EHA > 0)
# VN <- Qtd de EHC dado como EHC (qtd EHC < 0)
# FP <- Qtd de EHC dado como EHA (qtd EHC > 0)
# FN <- Qtd de EHA dado como EHC (qtd EHA < 0)
df_predict2_1 <- data.frame(df_predict2$.pred,rownames(df_predict2), categoria)
VP2 <- df_predict2_1[which(df_predict2_1$df_predict2..pred>0 & df_predict2_1$categoria=="EHA"),]
nrow(VP)
## [1] 91
VN2 <- df_predict1[which(df_predict2_1$df_predict2..pred<0 & df_predict2_1$categoria=="EHC"),]
nrow(VN)
## [1] 98
FP2 <- df_predict1[which(df_predict2_1$df_predict2..pred>0 & df_predict2_1$categoria=="EHC"),]
nrow(FP)
## [1] 0
FN2 <- df_predict1[which(df_predict2_1$df_predict2..pred<0 & df_predict2_1$categoria=="EHA"),]
nrow(FN)
## [1] 8
AUC2= 100*(nrow(VP2)+nrow(VN2))/(nrow(VP2)+nrow(VN2)+nrow(FP2)+nrow(FN2))
AUC2
## [1] 95.93909
Acurácia ficou a mesma.
df.marcadores3 <- data.frame(categoria,df[,c(28,32,61,14)])
df.marcadores3$categoria <- ifelse(df.marcadores3$categoria == "EHA", 1,
ifelse(df.marcadores3$categoria == "EHC", -1, df.marcadores3$categoria))
# Transformando a coluna 'categoria' de volta para numérico
df.marcadores3$categoria <- as.numeric(as.character(df.marcadores3$categoria))
# Verificando a transformação
str(df.marcadores3$categoria)
## num [1:197] 1 1 1 1 1 1 1 1 1 1 ...
#Regressão
lm_model3 <- linear_reg() %>%
set_engine('lm') %>%
set_mode('regression') %>%
fit(categoria ~ VOC25 * VOC29 * VOC59 * VOC11, data = df.marcadores3)
summary(lm_model3)
## Length Class Mode
## lvl 0 -none- NULL
## spec 7 linear_reg list
## fit 12 lm list
## preproc 1 -none- list
## elapsed 1 -none- list
## censor_probs 0 -none- list
#Predição
df_predict3 <- lm_model3 %>% predict(df)
rownames(df_predict3) <- rownames(df)
## Warning: Setting row names on a tibble is deprecated.
df_predict3
## # A tibble: 197 × 1
## .pred
## * <dbl>
## 1 1.10
## 2 0.563
## 3 0.957
## 4 0.804
## 5 1.01
## 6 0.898
## 7 1.33
## 8 0.537
## 9 0.881
## 10 1.13
## # ℹ 187 more rows
summary(lm_model3)
## Length Class Mode
## lvl 0 -none- NULL
## spec 7 linear_reg list
## fit 12 lm list
## preproc 1 -none- list
## elapsed 1 -none- list
## censor_probs 0 -none- list
#Plot
plot_ly(df, y = ~ lm_model3$fit$fitted.values, x = ~categoria,
colors = categoria,
type = 'scatter',
alpha = 0.2,
mode = 'markers',
name = 'Tips',
text = rownames(df_predict3)) %>%
layout(xaxis = list(title = "Categoria"),
yaxis = list(title = "Valor Predito"),
title = "PLS-DA")
# Cálculo da acurácia
# Sendo, o positivo a aditivação:
# VP <- Qtd de EHA dado como EHA (qtd EHA > 0)
# VN <- Qtd de EHC dado como EHC (qtd EHC < 0)
# FP <- Qtd de EHC dado como EHA (qtd EHC > 0)
# FN <- Qtd de EHA dado como EHC (qtd EHA < 0)
df_predict3_1 <- data.frame(df_predict3$.pred,rownames(df_predict3), categoria)
VP3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred>0 & df_predict3_1$categoria=="EHA"),]
nrow(VP3)
## [1] 91
VN3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred<0 & df_predict3_1$categoria=="EHC"),]
nrow(VN3)
## [1] 98
FP3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred>0 & df_predict3_1$categoria=="EHC"),]
nrow(FP3)
## [1] 0
FN3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred<0 & df_predict3_1$categoria=="EHA"),]
nrow(FN3)
## [1] 8
AUC3= 100*(nrow(VP3)+nrow(VN3))/(nrow(VP3)+nrow(VN3)+nrow(FP3)+nrow(FN3))
AUC3
## [1] 95.93909
Foi escolhido o Modelo 1 com três marcadores!